/*


!
!  Dalton, a molecular electronic structure program
!  Copyright (C) The Dalton Authors (see AUTHORS file for details).
!
!  This program is free software; you can redistribute it and/or
!  modify it under the terms of the GNU Lesser General Public
!  License version 2.1 as published by the Free Software Foundation.
!
!  This program is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
!  Lesser General Public License for more details.
!
!  If a copy of the GNU LGPL v2.1 was not distributed with this
!  code, you can obtain one at https://www.gnu.org/licenses/old-licenses/lgpl-2.1.en.html.
!

!

*/
/*-*-mode: C; c-indentation-style: "bsd"; c-basic-offset: 4; -*-*/
/* fun-pw91x2.c:

   REFERENCE: P.M.W. Gill, Mol. Phys., 89, 433 (1996).
   see also C. Adamo and V. Barone, J. Chem. Phys., 108, 664 (1998)
   and the original PW91 paper of J.P. Perdew, J.A. Chevary, S.H. Vosko, K.A. Jackson,
       M.R. Pederson, D.J. Singh and C. Fiolhais, Phys. Rev. B, 46, 6671 (1992). 

   PW91x exchange functional as parameterized by Gill. 
   Implemented by David Wilson (david.wilson@latrobe.edu.au), Jun 2005

   Automatically generated code implementing PW91X2 functional and
   its derivatives. It is generated by func-codegen.pl being a part of
   a "Automatic code generation framework for analytical functional
   derivative evaluation", Pawel Salek, 2005

    This functional is connected by making following changes:
    1. add "extern Functional pw91x2Functional;" to 'functionals.h'
    2. add "&pw91x2Functional," to 'functionals.c'
    3. add "fun-pw91x2.c" to 'Makefile.am', 'Makefile.in' or 'Makefile'.

    This functional has been generated from following input:
    ------ cut here -------
rho:  rhoa + rhob;
grad: sqrt(grada*grada + gradb*gradb + 2*gradab);
zeta: (rhoa-rhob)/(rhoa+rhob);

b:    0.0042;
c:    1.6455;
d:    4;
Ax:   -3/4*(6/%PI)^(1/3);
bet:  5*((36*%PI)^(-5/3));

xa:    grada/(rhoa^(4/3));
xb:    gradb/(rhob^(4/3));

F0a:   b*xa^2-(b-bet)*xa^2*exp(-c*xa^2)-10^(-6)*xa^4;
F0b:   b*xb^2-(b-bet)*xb^2*exp(-c*xb^2)-10^(-6)*xb^4;

F1a:   1+6*b*xa*asinh(xa)-(10^(-6)*xa^4)/Ax;
F1b:   1+6*b*xb*asinh(xb)-(10^(-6)*xb^4)/Ax;

Fa: F0a/F1a;
Fb: F0b/F1b;

Exa: rhoa^(4/3)*(Ax-Fa);
Exb: rhob^(4/3)*(Ax-Fb);

K(rhoa,rhob,grada,gradb,gradab):=(Exa+Exb);


    ------ cut here -------
*/

 
/* strictly conform to XOPEN ANSI C standard */
#if !defined(SYS_DEC)
/* XOPEN compliance is missing on old Tru64 4.0E Alphas and pow() prototype
 * is not specified. */
#define _XOPEN_SOURCE          500
#define _XOPEN_SOURCE_EXTENDED 1
#endif
#include <math.h>
#include <stddef.h>
#include "general.h"

#define __CVERSION__
 
#include "functionals.h"

/* INTERFACE PART */
static integer pw91x2_isgga(void) { return 1; } /* FIXME: detect! */
static integer pw91x2_read(const char *conf_line);
static real pw91x2_energy(const FunDensProp* dp);
static void pw91x2_first(FunFirstFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void pw91x2_second(FunSecondFuncDrv *ds, real factor,
                          const FunDensProp* dp);
static void pw91x2_third(FunThirdFuncDrv *ds,   real factor,
                         const FunDensProp* dp);
static void pw91x2_fourth(FunFourthFuncDrv *ds,   real factor,
                          const FunDensProp* dp);
 
Functional PW91x2Functional = {
  "PW91x2",       /* name */
  pw91x2_isgga,   /* gga-corrected */
   1,
  pw91x2_read,
  NULL,
  pw91x2_energy,
  pw91x2_first,
  pw91x2_second,
  pw91x2_third,
  pw91x2_fourth
};
 
/* IMPLEMENTATION PART */
static integer
pw91x2_read(const char *conf_line)
{
    fun_set_hf_weight(0);
    return 1;
}

static real
pw91x2_energy(const FunDensProp *dp)
{
    real res;
    real rhoa = dp->rhoa, rhob = dp->rhob;
    real grada = dp->grada, gradb = dp->gradb, gradab = dp->gradab;

    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17;

    t1 = pow(rhoa,1.333333333333333);
    t2 = pow(6.0,0.333333333333333);
    t3 = pow(3.141592653589793,0.333333333333333);
    t4 = -0.75*t2/t3;
    t5 = 1/t2;
    t6 = pow(grada,4.0);
    t7 = 1/pow(rhoa,5.333333333333333);
    t8 = 1/t1;
    t9 = pow(grada,2.0);
    t10 = 1/pow(rhoa,2.666666666666667);
    t11 = 0.0042-0.138888888888889/(pow(3.141592653589793,
        1.666666666666667)*pow(36.0,0.666666666666667));
    t12 = pow(rhob,1.333333333333333);
    t13 = pow(gradb,4.0);
    t14 = 1/pow(rhob,5.333333333333333);
    t15 = 1/t12;
    t16 = pow(gradb,2.0);
    t17 = 1/pow(rhob,2.666666666666667);

   /* code */
    res = t1*(t4-1.0*(-1.0*t10*t11*t9/pow(2.718281828459045,
        1.6455*t9*t10)-1.0E-6*t6*t7+0.0042*t9*t10)/(0.0252*grada*asinh(grada*
        t8)*t8+1.3333333333333334E-6*t3*t5*t6*t7+1.0))+t12*(t4-1.0*
        (-1.0*t11*t16*t17/pow(2.718281828459045,1.6455*t16*t17)+0.0042*
        t16*t17-1.0E-6*t13*t14)/(1.3333333333333334E-6*t13*t14*t3*
        t5+0.0252*gradb*asinh(gradb*t15)*t15+1.0));

    return res;
}

static void
pw91x2_first_helper(real rhoa, real grada, real *res)
{    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20;

    t1 = pow(rhoa,1.333333333333333);
    t2 = pow(6.0,0.333333333333333);
    t3 = 1/t2;
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/t1;
    t8 = asinh(grada*t7);
    t9 = 0.0252*grada*t8*t7+1.3333333333333334E-6*t3*t4*t5*
        t6+1.0;
    t10 = 1/t9;
    t11 = 1/pow(rhoa,6.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,3.666666666666667);
    t14 = 0.0042-0.138888888888889/(pow(3.141592653589793,
        1.666666666666667)*pow(36.0,0.666666666666667));
    t15 = 1/pow(rhoa,2.666666666666667);
    t16 = 1/pow(2.718281828459045,1.6455*t12*t15);
    t17 = 1/sqrt(t12*t15+1.0);
    t18 = 1/pow(t9,2.0);
    t19 = -1.0E-6*t5*t6-1.0*t12*t14*t15*t16+0.0042*t12*t15;
    t20 = pow(grada,
        3.0);

   /* code */
    res[0] = t1*(t18*t19*(-0.0336*t8*grada/pow(rhoa,2.333333333333333)-
        7.111111111111111E-6*t11*t3*t4*t5-0.0336*t12*t17*t13)-1.0*
        t10*(5.333333333333334E-6*t11*t5+2.666666666666667*t12*t13*
        t14*t16-4.388*t14*t5*t11*t16-0.0112*t12*t13))+1.333333333333333*
        (-0.75*t2/t4-1.0*t10*t19)*pow(rhoa,0.333333333333333);
    res[1] = t1*(t18*t19*(0.0252*t8*t7+5.333333333333334E-6*
        t20*t3*t4*t6+0.0252*grada*t17*t15)-1.0*t10*(-2.0*t14*t15*t16*
        grada-4.0E-6*t20*t6+3.291*t14*t20*t6*t16+0.0084*grada*t15));
}

static void
pw91x2_first(FunFirstFuncDrv *ds, real factor, const FunDensProp *dp)
{
    real res[2];

    pw91x2_first_helper(dp->rhoa, dp->grada, res);
   /* Final assignment */
    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x2_first_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

}

static void
pw91x2_second_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;

    t1 = pow(rhoa,1.333333333333333);
    t2 = pow(6.0,0.333333333333333);
    t3 = 1/t2;
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/t1;
    t8 = asinh(grada*t7);
    t9 = 0.0252*grada*t8*t7+1.3333333333333334E-6*t3*t4*t5*
        t6+1.0;
    t10 = 1/t9;
    t11 = 1/pow(rhoa,6.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,3.666666666666667);
    t14 = 0.0042-0.138888888888889/(pow(3.141592653589793,
        1.666666666666667)*pow(36.0,0.666666666666667));
    t15 = 1/pow(rhoa,2.666666666666667);
    t16 = 1/pow(2.718281828459045,1.6455*t12*t15);
    t17 = 5.333333333333334E-6*t11*t5+2.666666666666667*t12*
        t13*t14*t16-4.388*t14*t5*t11*t16-0.0112*t12*t13;
    t18 = sqrt(t12*t15+1.0);
    t19 = 1/t18;
    t20 = 1/pow(rhoa,2.333333333333333);
    t21 = -7.111111111111111E-6*t11*t3*t4*t5-0.0336*grada*
        t8*t20-0.0336*t12*t19*t13;
    t22 = 1/pow(t9,2.0);
    t23 = -1.0E-6*t5*t6-1.0*t12*t14*t15*t16+0.0042*t12*t15;
    t24 = t21*
        t22*t23-1.0*t10*t17;
    t25 = pow(rhoa,0.333333333333333);
    t26 = -0.75*t2/t4-1.0*t10*t23;
    t27 = pow(grada,3.0);
    t28 = -2.0*t14*t15*t16*grada-4.0E-6*t27*t6+3.291*t14*
        t27*t6*t16+0.0084*grada*t15;
    t29 = 0.0252*t8*t7+5.333333333333334E-6*t27*t3*t4*t6+
        0.0252*grada*t19*t15;
    t30 = t29*t22*t23-1.0*t10*t28;
    t31 = 1/pow(rhoa,7.333333333333333);
    t32 = 1/pow(rhoa,4.666666666666667);
    t33 = 1/pow(t9,3.0);
    t34 = 1/pow(t18,3.0);

   /* code */
    res[0] = 1.333333333333333*t25*t26+t1*t24;
    res[1] = t1*t30;
    res[2] = t1*(-1.0*t10*(-19.254544*t14*t16*pow(grada,6.0)/
        pow(rhoa,10.0)-3.377777777777778E-5*t31*t5-9.777777777777779*
        t12*t14*t16*t32+0.041066666666667*t12*t32+39.492*t14*t5*t31*
        t16)+t22*t23*(0.0784*t8*grada/pow(rhoa,3.333333333333333)+
        4.503703703703703E-5*t3*t31*t4*t5+0.168*t12*t19*t32-0.0448*
        t5*t34*t31)-2.0*pow(t21,2.0)*t23*t33+2.0*t17*t21*t22)+0.444444444444444*
        t26/pow(rhoa,0.666666666666667)+2.666666666666667*t24*t25;
    res[3] = t1*(-1.0*t10*(14.440908*t14*t16*pow(grada,5.0)/
        pow(rhoa,9.0)+5.333333333333333*t13*t14*t16*grada+2.1333333333333335E-5*
        t11*t27-26.328*t14*t27*t11*t16-0.0224*grada*t13)+t22*t23*(-
        2.8444444444444444E-5*t11*t27*t3*t4-0.0336*t8*t20-0.1008*grada*
        t19*t13+0.0336*t27*t34*t11)-2.0*t21*t23*t29*t33+t21*t22*t28+
        t29*t22*t17)+1.333333333333333*t25*t30;
    res[4] = t1*(-1.0*t10*(-10.830681*t14*t16*t5/pow(rhoa,
        8.0)-1.2E-5*t12*t6+16.455*t14*t12*t6*t16-2.0*t14*t15*t16+0.0084*
        t15)+t22*t23*(1.6E-5*t12*t3*t4*t6-0.0252*t12*t34*t6+0.0504*
        t19*t15)-2.0*t23*pow(t29,2.0)*t33+2.0*t22*t28*t29);

}

static void
pw91x2_second(FunSecondFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[5];
 
    pw91x2_second_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x2_second_helper(dp->rhob, dp->gradb, res);
    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

}

static void
pw91x2_third_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57;

    t1 = pow(rhoa,1.333333333333333);
    t2 = pow(6.0,0.333333333333333);
    t3 = 1/t2;
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/t1;
    t8 = asinh(grada*t7);
    t9 = 0.0252*grada*t8*t7+1.3333333333333334E-6*t3*t4*t5*
        t6+1.0;
    t10 = 1/t9;
    t11 = 1/pow(rhoa,6.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,3.666666666666667);
    t14 = 0.0042-0.138888888888889/(pow(3.141592653589793,
        1.666666666666667)*pow(36.0,0.666666666666667));
    t15 = 1/pow(rhoa,2.666666666666667);
    t16 = 1/pow(2.718281828459045,1.6455*t12*t15);
    t17 = 5.333333333333334E-6*t11*t5+2.666666666666667*t12*
        t13*t14*t16-4.388*t14*t5*t11*t16-0.0112*t12*t13;
    t18 = sqrt(t12*t15+1.0);
    t19 = 1/t18;
    t20 = 1/pow(rhoa,2.333333333333333);
    t21 = -7.111111111111111E-6*t11*t3*t4*t5-0.0336*grada*
        t8*t20-0.0336*t12*t19*t13;
    t22 = 1/pow(t9,2.0);
    t23 = -1.0E-6*t5*t6-1.0*t12*t14*t15*t16+0.0042*t12*t15;
    t24 = t21*
        t22*t23-1.0*t10*t17;
    t25 = pow(rhoa,0.333333333333333);
    t26 = -0.75*t2/t4-1.0*t10*t23;
    t27 = pow(grada,3.0);
    t28 = -2.0*t14*t15*t16*grada-4.0E-6*t27*t6+3.291*t14*
        t27*t6*t16+0.0084*grada*t15;
    t29 = 0.0252*t8*t7+5.333333333333334E-6*t27*t3*t4*t6+
        0.0252*grada*t19*t15;
    t30 = t29*t22*t23-1.0*t10*t28;
    t31 = 1/pow(rhoa,7.333333333333333);
    t32 = 1/pow(rhoa,4.666666666666667);
    t33 = pow(grada,6.0);
    t34 = 1/pow(rhoa,10.0);
    t35 = -3.377777777777778E-5*t31*t5-9.777777777777779*
        t12*t14*t16*t32+0.041066666666667*t12*t32-19.254544*t14*t33*
        t34*t16+39.492*t14*t5*t31*t16;
    t36 = pow(t21,2.0);
    t37 = 1/pow(t9,3.0);
    t38 = 1/pow(t18,3.0);
    t39 = 1/pow(rhoa,3.333333333333333);
    t40 = 4.503703703703703E-5*t3*t31*t4*t5+0.0784*grada*
        t8*t39+0.168*t12*t19*t32-0.0448*t5*t38*t31;
    t41 = -2.0*t23*t36*t37-1.0*t10*t35+t40*t22*t23+2.0*t17*
        t21*t22;
    t42 = 1/pow(rhoa,0.666666666666667);
    t43 = pow(grada,5.0);
    t44 = 1/pow(rhoa,9.0);
    t45 = 5.333333333333333*t13*t14*t16*grada+2.1333333333333335E-5*
        t11*t27+14.440908*t14*t43*t44*t16-26.328*t14*t27*t11*t16-0.0224*
        grada*t13;
    t46 = -2.8444444444444444E-5*t11*t27*t3*t4-0.0336*t8*
        t20-0.1008*grada*t19*t13+0.0336*t27*t38*t11;
    t47 = -1.0*t10*t45-2.0*t21*t23*t29*t37+t21*t22*t28+t46*
        t22*t23+t29*t22*t17;
    t48 = 1/pow(rhoa,8.0);
    t49 = -1.2E-5*t12*t6+16.455*t14*t12*t6*t16-10.830681*
        t14*t5*t48*t16-2.0*t14*t15*t16+0.0084*t15;
    t50 = pow(t29,2.0);
    t51 = 1.6E-5*t12*t3*t4*t6-0.0252*t12*t38*t6+0.0504*t19*
        t15;
    t52 = -2.0*t23*t37*t50-1.0*t10*t49+2.0*t22*t28*t29+t51*
        t22*t23;
    t53 = 1/pow(rhoa,8.333333333333334);
    t54 = 1/pow(rhoa,5.666666666666667);
    t55 = 1/pow(rhoa,11.0);
    t56 = 1/pow(t9,4.0);
    t57 = 1/pow(t18,5.0);

   /* code */
    res[0] = 1.333333333333333*t25*t26+t1*t24;
    res[1] = t1*t30;
    res[2] = 0.444444444444444*t26*t42+t1*t41+2.666666666666667*
        t24*t25;
    res[3] = t1*t47+1.333333333333333*t25*t30;
    res[4] = t1*t52;
    res[5] = t1*(-1.0*t10*(-84.488939072*t14*t16*pow(grada,
        8.0)/pow(rhoa,13.66666666666667)+45.62962962962963*t12*t14*
        t16*t54-0.191644444444444*t12*t54+2.477037037037037E-4*t5*
        t53+365.836336*t14*t33*t55*t16-332.5128888888888*t14*t5*t53*
        t16)+t22*t23*(-0.261333333333333*t8*grada/pow(rhoa,4.333333333333333)-
        0.1792*t33*t57*t55-0.888533333333333*t12*t19*t54-3.302716049382716E-4*
        t3*t4*t5*t53+0.552533333333333*t5*t38*t53)+6.0*pow(t21,3.0)*
        t23*t56-6.0*t21*t23*t37*t40+3.0*t17*t22*t40-6.0*t17*t36*t37+
        3.0*t21*t22*t35)-0.296296296296296*t26/pow(rhoa,1.666666666666667)+
        1.333333333333333*t24*t42+4.0*t25*t41;
    res[6] = t1*(-1.0*t10*(63.366704304*t14*t16*pow(grada,
        7.0)/pow(rhoa,12.66666666666667)-19.55555555555556*t14*t16*
        t32*grada+0.082133333333333*grada*t32-1.351111111111111E-4*
        t27*t31-245.495436*t14*t43*t34*t16+190.1466666666666*t14*t27*
        t31*t16)+6.0*t23*t29*t36*t56-4.0*t21*t23*t37*t46+2.0*t17*t22*
        t46+2.0*t21*t22*t45-2.0*t23*t29*t37*t40+t22*t23*(1.8014814814814814E-4*
        t27*t3*t31*t4+0.0784*t8*t39+0.1344*t43*t57*t34+0.4144*grada*
        t19*t32-0.3472*t27*t38*t31)-2.0*t28*t36*t37-4.0*t17*t21*t29*
        t37+t29*t22*t35+t40*t22*t28)+2.666666666666667*t25*t47+0.444444444444444*
        t30*t42;
    res[7] = t1*(-1.0*t10*(-47.525028228*t14*t16*t33/pow(rhoa,
        11.66666666666667)+158.849988*t14*t5*t44*t16+5.333333333333333*
        t13*t14*t16-96.536*t14*t12*t11*t16-0.0224*t13+6.4E-5*t11*t12)+
        6.0*t21*t23*t50*t56-2.0*t21*t23*t37*t51-2.0*t17*t37*t50+t21*
        t22*t49-4.0*t23*t29*t37*t46+2.0*t22*t28*t46+2.0*t22*t29*t45+
        t22*t23*(-0.1008*t5*t57*t44-8.533333333333334E-5*t11*t12*t3*
        t4-0.1344*t19*t13+0.2016*t12*t38*t11)-4.0*t21*t28*t29*t37+
        t51*t22*t17)+1.333333333333333*t25*t52;
    res[8] = t1*(-1.0*t10*(35.643771171*t14*t16*t43/pow(rhoa,
        10.66666666666667)-2.4E-5*t6*grada+39.492*t14*grada*t6*t16-
        97.47612899999999*t14*t27*t48*t16)+t22*t23*(3.2E-5*t3*t4*t6*
        grada-0.1008*grada*t38*t6+0.0756*t27*t57*t48)+6.0*t23*pow(t29,
        3.0)*t56-6.0*t23*t29*t37*t51+3.0*t22*t28*t51-6.0*t28*t37*t50+
        3.0*t22*t29*t49);

}

static void
pw91x2_third(FunThirdFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[9];
 
    pw91x2_third_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x2_third_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

}

static void
pw91x2_fourth_helper(real rhoa, real grada, real *res)
{
    real t1, t2, t3, t4, t5, t6, t7, t8, t9, t10;
    real t11, t12, t13, t14, t15, t16, t17, t18;
    real t19, t20, t21, t22, t23, t24, t25, t26;
    real t27, t28, t29, t30, t31, t32, t33, t34;
    real t35, t36, t37, t38, t39, t40, t41, t42;
    real t43, t44, t45, t46, t47, t48, t49, t50;
    real t51, t52, t53, t54, t55, t56, t57, t58;
    real t59, t60, t61, t62, t63, t64, t65, t66;
    real t67, t68, t69, t70, t71, t72, t73, t74;
    real t75, t76, t77, t78, t79, t80, t81, t82;
    real t83, t84, t85;

    t1 = pow(rhoa,1.333333333333333);
    t2 = pow(6.0,0.333333333333333);
    t3 = 1/t2;
    t4 = pow(3.141592653589793,0.333333333333333);
    t5 = pow(grada,4.0);
    t6 = 1/pow(rhoa,5.333333333333333);
    t7 = 1/t1;
    t8 = asinh(grada*t7);
    t9 = 0.0252*grada*t8*t7+1.3333333333333334E-6*t3*t4*t5*
        t6+1.0;
    t10 = 1/t9;
    t11 = 1/pow(rhoa,6.333333333333333);
    t12 = pow(grada,2.0);
    t13 = 1/pow(rhoa,3.666666666666667);
    t14 = 0.0042-0.138888888888889/(pow(3.141592653589793,
        1.666666666666667)*pow(36.0,0.666666666666667));
    t15 = 1/pow(rhoa,2.666666666666667);
    t16 = 1/pow(2.718281828459045,1.6455*t12*t15);
    t17 = 5.333333333333334E-6*t11*t5+2.666666666666667*t12*
        t13*t14*t16-4.388*t14*t5*t11*t16-0.0112*t12*t13;
    t18 = sqrt(t12*t15+1.0);
    t19 = 1/t18;
    t20 = 1/pow(rhoa,2.333333333333333);
    t21 = -7.111111111111111E-6*t11*t3*t4*t5-0.0336*grada*
        t8*t20-0.0336*t12*t19*t13;
    t22 = 1/pow(t9,2.0);
    t23 = -1.0E-6*t5*t6-1.0*t12*t14*t15*t16+0.0042*t12*t15;
    t24 = t21*
        t22*t23-1.0*t10*t17;
    t25 = pow(rhoa,0.333333333333333);
    t26 = -0.75*t2/t4-1.0*t10*t23;
    t27 = pow(grada,3.0);
    t28 = -2.0*t14*t15*t16*grada-4.0E-6*t27*t6+3.291*t14*
        t27*t6*t16+0.0084*grada*t15;
    t29 = 0.0252*t8*t7+5.333333333333334E-6*t27*t3*t4*t6+
        0.0252*grada*t19*t15;
    t30 = t29*t22*t23-1.0*t10*t28;
    t31 = 1/pow(rhoa,7.333333333333333);
    t32 = 1/pow(rhoa,4.666666666666667);
    t33 = pow(grada,6.0);
    t34 = 1/pow(rhoa,10.0);
    t35 = -3.377777777777778E-5*t31*t5-9.777777777777779*
        t12*t14*t16*t32+0.041066666666667*t12*t32-19.254544*t14*t33*
        t34*t16+39.492*t14*t5*t31*t16;
    t36 = pow(t21,2.0);
    t37 = 1/pow(t9,3.0);
    t38 = 1/pow(t18,3.0);
    t39 = 1/pow(rhoa,3.333333333333333);
    t40 = 4.503703703703703E-5*t3*t31*t4*t5+0.0784*grada*
        t8*t39+0.168*t12*t19*t32-0.0448*t5*t38*t31;
    t41 = -2.0*t23*t36*t37-1.0*t10*t35+t40*t22*t23+2.0*t17*
        t21*t22;
    t42 = 1/pow(rhoa,0.666666666666667);
    t43 = pow(grada,5.0);
    t44 = 1/pow(rhoa,9.0);
    t45 = 5.333333333333333*t13*t14*t16*grada+2.1333333333333335E-5*
        t11*t27+14.440908*t14*t43*t44*t16-26.328*t14*t27*t11*t16-0.0224*
        grada*t13;
    t46 = -2.8444444444444444E-5*t11*t27*t3*t4-0.0336*t8*
        t20-0.1008*grada*t19*t13+0.0336*t27*t38*t11;
    t47 = -1.0*t10*t45-2.0*t21*t23*t29*t37+t21*t22*t28+t46*
        t22*t23+t29*t22*t17;
    t48 = 1/pow(rhoa,8.0);
    t49 = -1.2E-5*t12*t6+16.455*t14*t12*t6*t16-10.830681*
        t14*t5*t48*t16-2.0*t14*t15*t16+0.0084*t15;
    t50 = pow(t29,2.0);
    t51 = 1.6E-5*t12*t3*t4*t6-0.0252*t12*t38*t6+0.0504*t19*
        t15;
    t52 = -2.0*t23*t37*t50-1.0*t10*t49+2.0*t22*t28*t29+t51*
        t22*t23;
    t53 = 1/pow(rhoa,8.333333333333334);
    t54 = 1/pow(rhoa,5.666666666666667);
    t55 = pow(grada,8.0);
    t56 = 1/pow(rhoa,13.66666666666667);
    t57 = 1/pow(rhoa,11.0);
    t58 = 45.62962962962963*t12*t14*t16*t54-0.191644444444444*
        t12*t54+2.477037037037037E-4*t5*t53+365.836336*t14*t33*t57*
        t16-84.488939072*t14*t55*t56*t16-332.5128888888888*t14*t5*
        t53*t16;
    t59 = pow(t21,3.0);
    t60 = 1/pow(t9,4.0);
    t61 = 1/pow(t18,5.0);
    t62 = 1/pow(rhoa,4.333333333333333);
    t63 = -0.261333333333333*grada*t8*t62-0.1792*t33*t61*
        t57-0.888533333333333*t12*t19*t54-3.302716049382716E-4*t3*
        t4*t5*t53+0.552533333333333*t5*t38*t53;
    t64 = 6.0*t23*t59*t60-1.0*t10*t58-6.0*t21*t23*t37*t40+
        3.0*t17*t22*t40-6.0*t17*t36*t37+3.0*t21*t22*t35+t63*t22*t23;
    t65 = 1/
        pow(rhoa,1.666666666666667);
    t66 = pow(grada,7.0);
    t67 = 1/pow(rhoa,12.66666666666667);
    t68 = -19.55555555555556*t14*t16*t32*grada+0.082133333333333*
        grada*t32-1.351111111111111E-4*t27*t31+63.366704304*t14*t66*
        t67*t16-245.495436*t14*t43*t34*t16+190.1466666666666*t14*t27*
        t31*t16;
    t69 = 1.8014814814814814E-4*t27*t3*t31*t4+0.0784*t8*t39+
        0.1344*t43*t61*t34+0.4144*grada*t19*t32-0.3472*t27*t38*t31;
    t70 = -
        1.0*t10*t68+6.0*t23*t29*t36*t60-4.0*t21*t23*t37*t46+2.0*t17*
        t22*t46+2.0*t21*t22*t45-2.0*t23*t29*t37*t40-2.0*t28*t36*t37-
        4.0*t17*t21*t29*t37+t29*t22*t35+t40*t22*t28+t69*t22*t23;
    t71 = 1/
        pow(rhoa,11.66666666666667);
    t72 = -47.525028228*t14*t33*t71*t16+158.849988*t14*t5*
        t44*t16+5.333333333333333*t13*t14*t16-96.536*t14*t12*t11*t16-
        0.0224*t13+6.4E-5*t11*t12;
    t73 = -0.1008*t5*t61*t44-8.533333333333334E-5*t11*t12*
        t3*t4-0.1344*t19*t13+0.2016*t12*t38*t11;
    t74 = -1.0*t10*t72+6.0*t21*t23*t50*t60-2.0*t21*t23*t37*
        t51-2.0*t17*t37*t50+t21*t22*t49-4.0*t23*t29*t37*t46+2.0*t22*
        t28*t46+2.0*t22*t29*t45-4.0*t21*t28*t29*t37+t73*t22*t23+t51*
        t22*t17;
    t75 = 1/pow(rhoa,10.66666666666667);
    t76 = -2.4E-5*t6*grada+35.643771171*t14*t43*t75*t16+39.492*
        t14*grada*t6*t16-97.47612899999999*t14*t27*t48*t16;
    t77 = pow(t29,3.0);
    t78 = 3.2E-5*t3*t4*t6*grada-0.1008*grada*t38*t6+0.0756*
        t27*t61*t48;
    t79 = 6.0*t23*t60*t77-1.0*t10*t76-6.0*t23*t29*t37*t51+
        3.0*t22*t28*t51-6.0*t28*t37*t50+3.0*t22*t29*t49+t78*t22*t23;
    t80 = 1/
        pow(rhoa,9.333333333333334);
    t81 = 1/pow(rhoa,6.666666666666667);
    t82 = 1/pow(rhoa,14.66666666666667);
    t83 = 1/pow(rhoa,12.0);
    t84 = 1/pow(t9,5.0);
    t85 = 1/pow(t18,7.0);

   /* code */
    res[0] = 1.333333333333333*t25*t26+t1*t24;
    res[1] = t1*t30;
    res[2] = 0.444444444444444*t26*t42+t1*t41+2.666666666666667*
        t24*t25;
    res[3] = t1*t47+1.333333333333333*t25*t30;
    res[4] = t1*t52;
    res[5] = -0.296296296296296*t26*t65+t1*t64+1.333333333333333*
        t24*t42+4.0*t25*t41;
    res[6] = t1*t70+2.666666666666667*t25*t47+0.444444444444444*
        t30*t42;
    res[7] = t1*t74+1.333333333333333*t25*t52;
    res[8] = t1*t79;
    res[9] = t1*(-1.0*t10*(-370.737464647936*t14*t16*pow(grada,
        10.0)/pow(rhoa,17.33333333333333)-258.5679012345679*t12*t14*
        t16*t81+1.085985185185185*t12*t81-0.002064197530864*t5*t80-
        5483.266252444444*t14*t33*t83*t16+2759.972009685333*t14*t55*
        t82*t16+2971.163555555555*t14*t5*t80*t16)-24.0*pow(t21,4.0)*
        t23*t84+t22*t23*(4.181333333333333*t33*t61*t83-1.194666666666667*
        t55*t85*t82+5.383466666666667*t12*t19*t81+0.002752263374486*
        t3*t4*t5*t80-5.789155555555556*t5*t38*t80+1.132444444444445*
        grada*t8*t6)-8.0*t21*t23*t37*t63+4.0*t17*t22*t63+24.0*t17*
        t59*t60+36.0*t23*t36*t40*t60+4.0*t21*t22*t58-6.0*t23*t37*pow(t40,
        2.0)-24.0*t17*t21*t37*t40+6.0*t22*t35*t40-12.0*t35*t36*t37)-
        1.185185185185185*t24*t65+5.333333333333333*t25*t64+2.666666666666667*
        t41*t42+0.493827160493827*t15*t26;
    res[10] = t1*(-1.0*t10*(278.053098485952*t14*t16*pow(grada,
        9.0)/pow(rhoa,16.33333333333333)+91.25925925925925*t14*t16*
        t54*grada-0.383288888888889*grada*t54+9.908148148148148E-4*
        t27*t53+3289.317933333333*t14*t43*t57*t16-1879.878894352*t14*
        t66*t56*t16-1480.218666666666*t14*t27*t53*t16)-24.0*t23*t29*
        t59*t84-6.0*t21*t23*t37*t69+3.0*t17*t22*t69+3.0*t21*t22*t68-
        2.0*t23*t29*t37*t63+t22*t23*(-0.261333333333333*t8*t62-2.7328*
        t43*t61*t57+0.896*t66*t85*t56-2.0384*grada*t19*t54-0.001321086419753*
        t27*t3*t4*t53+3.098666666666666*t27*t38*t53)+6.0*t28*t59*t60+
        18.0*t23*t36*t46*t60+18.0*t21*t23*t29*t40*t60+18.0*t17*t29*
        t36*t60+t29*t22*t58-6.0*t23*t37*t40*t46-12.0*t17*t21*t37*t46+
        3.0*t22*t35*t46+3.0*t22*t40*t45-6.0*t36*t37*t45-6.0*t17*t29*
        t37*t40-6.0*t21*t28*t37*t40-6.0*t21*t29*t35*t37+t63*t22*t28)+
        4.0*t25*t70-0.296296296296296*t30*t65+1.333333333333333*t42*
        t47;
    res[11] = t1*(-1.0*t10*(-208.539823864464*t14*t16*t55/
        pow(rhoa,15.33333333333333)-19.55555555555556*t14*t16*t32+
        0.082133333333333*t32-4.053333333333333E-4*t12*t31+1251.492410004*
        t14*t33*t67*t16-1853.24986*t14*t5*t34*t16+634.7973333333333*
        t14*t12*t31*t16)-24.0*t23*t36*t50*t84-4.0*t21*t23*t37*t73+
        2.0*t17*t22*t73+2.0*t21*t22*t72-4.0*t23*t29*t37*t69+2.0*t22*
        t28*t69+2.0*t22*t29*t68+t22*t23*(-0.672*t33*t85*t67+5.404444444444444E-4*
        t12*t3*t31*t4+1.7136*t5*t61*t34+0.4928*t19*t32-1.456*t12*t38*
        t31)+6.0*t23*t36*t51*t60+6.0*t23*t40*t50*t60+12.0*t17*t21*
        t50*t60+24.0*t21*t23*t29*t46*t60+12.0*t28*t29*t36*t60-2.0*
        t23*t37*t40*t51-4.0*t17*t21*t37*t51-2.0*t35*t37*t50-2.0*t36*
        t37*t49+t40*t22*t49-4.0*t23*t37*pow(t46,2.0)+4.0*t22*t45*t46-
        8.0*t17*t29*t37*t46-8.0*t21*t28*t37*t46-8.0*t21*t29*t37*t45-
        4.0*t28*t29*t37*t40+t51*t22*t35)+2.666666666666667*t25*t74+
        0.444444444444444*t42*t52;
    res[12] = t1*(-1.0*t10*(156.404867898348*t14*t16*t66/
        pow(rhoa,14.33333333333333)+1.28E-4*t11*grada-807.9254798759999*
        t14*t43*t71*t16+953.099928*t14*t27*t44*t16-210.624*t14*grada*
        t11*t16)+t22*t23*(-1.7066666666666668E-4*t11*t3*t4*grada+0.504*
        t43*t85*t71-1.008*t27*t61*t44+0.5376*grada*t38*t11)-24.0*t21*
        t23*t77*t84-2.0*t21*t23*t37*t78+6.0*t17*t60*t77+t21*t22*t76-
        6.0*t23*t29*t37*t73+3.0*t22*t28*t73+3.0*t22*t29*t72+18.0*t21*
        t23*t29*t51*t60+18.0*t23*t46*t50*t60+18.0*t21*t28*t50*t60-
        6.0*t23*t37*t46*t51+3.0*t22*t45*t51-6.0*t17*t29*t37*t51-6.0*
        t21*t28*t37*t51-6.0*t37*t45*t50+3.0*t22*t46*t49-6.0*t21*t29*
        t37*t49-12.0*t28*t29*t37*t46+t78*t22*t17)+1.333333333333333*
        t25*t79;
    res[13] = t1*(-1.0*t10*(-117.303650923761*t14*t16*t33/
        pow(rhoa,13.33333333333333)-2.4E-5*t6+499.0127963939999*t14*
        t5*t75*t16+39.492*t14*t6*t16-422.3965589999999*t14*t12*t48*
        t16)-24.0*t23*pow(t29,4.0)*t84-8.0*t23*t29*t37*t78+4.0*t22*
        t28*t78+24.0*t28*t60*t77+4.0*t22*t29*t76+t22*t23*(-0.378*t5*
        t85*t75+3.2E-5*t3*t4*t6-0.1008*t38*t6+0.5292*t12*t61*t48)+
        36.0*t23*t50*t51*t60-6.0*t23*t37*pow(t51,2.0)+6.0*t22*t49*
        t51-24.0*t28*t29*t37*t51-12.0*t37*t49*t50);

}

static void
pw91x2_fourth(FunFourthFuncDrv *ds, real factor, const FunDensProp* dp)
{
    real res[14];
 
    pw91x2_fourth_helper(dp->rhoa, dp->grada, res);

    ds->df1000 += factor*res[0];
    ds->df0010 += factor*res[1];

    ds->df2000 += factor*res[2];
    ds->df1010 += factor*res[3];
    ds->df0020 += factor*res[4];

    ds->df3000 += factor*res[5];
    ds->df2010 += factor*res[6];
    ds->df1020 += factor*res[7];
    ds->df0030 += factor*res[8];

    ds->df4000 += factor*res[9];
    ds->df3010 += factor*res[10];
    ds->df2020 += factor*res[11];
    ds->df1030 += factor*res[12];
    ds->df0040 += factor*res[13];


    if(fabs(dp->rhoa-dp->rhob)>1e-13 ||
       fabs(dp->grada-dp->gradb)>1e-13)
        pw91x2_fourth_helper(dp->rhob, dp->gradb, res);

    ds->df0100 += factor*res[0];
    ds->df0001 += factor*res[1];

    ds->df0200 += factor*res[2];
    ds->df0101 += factor*res[3];
    ds->df0002 += factor*res[4];

    ds->df0300 += factor*res[5];
    ds->df0201 += factor*res[6];
    ds->df0102 += factor*res[7];
    ds->df0003 += factor*res[8];

    ds->df0400 += factor*res[9];
    ds->df0301 += factor*res[10];
    ds->df0202 += factor*res[11];
    ds->df0103 += factor*res[12];
    ds->df0004 += factor*res[13];

}
